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airfoil chord, m 

skin friction coefficient 

specific heat, J/kg 

fraction of surface which is wetted. 

laminar convective heat transfer coefficient, W/m 2o K 

turbulent convective heat transfer coefficient, W/m 2o K 

equivalent sand grain roughness, m 

liquid water content, g/m 

median volume droplet diameter, pm 

static pressure, Pa 

pressure at the edge of the boundary layer, Pa 

turbulent Prantl No. 

relative humidity, % 

local roughness Reynolds No. 

surface distance from the stagnation point, m 

roughness Stanton No. 

temperature at the edge of the boundary layer, °C 

static air temperature, °C 

velocity at the edge of the boundary layer, m/s 

velocity at roughness height k, m/s 

shear velocity, m/s 

freestream velocity, m/s 

angle of attack, degrees 

momentum thickness, m 
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V 

X 

a 

x 


viscosity, Ns/m 2 

thermal conductivity of air, W/m°K 
viscosity, kg/m/s 
surface tension, N/m 
shear stress, N/m 2 


2.0 SUMMARY 


This report presents results from a joint study by DRA, NASA, and ONERA for the purpose of comparing, 
improving, and validating the aircraft icing computer codes developed by each agency. These codes are of 
three kinds: 1) water droplet trajectory prediction, 2) ice accretion modeling, and 3) transient electrothermal 
deicer analysis. In this joint study, the agencies compared their code predictions with each other and with 
experimental results. These comparison exercises were published in three technical reports, each with joint 
authorship. DRA published and had first authorship of Part I - Droplet Trajectory Calculations, NASA of Part II - 
Ice Accretion Prediction, and ONERA of Part III - Electrothermal Deicer Analysis. The results cover work done 
during the period from August 1 986 to late 1 991 . As a result, all of the information in this report is dated. Where 
necessary, current information is provided to show the direction of current research. 

In this present report on ice accretion, each agency predicted ice shapes on two dimensional airfoils under 
icing conditions for which experimental ice shapes were available. In general, all three codes did a reasonable 
job of predicting the measured ice shapes. For any given experimental condition, one of the three codes pre- 
dicted the general ice features (i.e., shape, impingement limits, mass of ice) somewhat better than did the 
other two. However, no single code consistently did better than the other two over the full range of conditions 
examined, which included rime, mixed, and glaze ice conditions. In several of the cases, DRA showed that the 
user’s knowledge of icing can significantly improve the accuracy of the code prediction. 

Rime ice predictions were reasonably accurate and consistent among the codes, because droplets freeze 
on impact and the freezing model is simple. Glaze ice predictions were less accurate and less consistent 
among the codes, because the freezing model is more complex and is critically dependent upon unsubstanti- 
ated heat transfer and surface roughness models. Thus, heat transfer prediction methods used in the codes 
became the subject for a separate study in this report to compare predicted heat transfer coefficients with a 
limited experimental database of heat transfer coefficients for cylinders with simulated glaze and rime ice 
shapes. The codes did a good job of predicting heat transfer coefficients near the stagnation region of the ice 
shapes. But in the region of the ice horns, all three codes predicted heat transfer coefficients considerably 
higher than the measured values. 

An important conclusion of this study is that further research is needed to understand the finer details of 
the glaze ice accretion process and to develop improved glaze ice accretion models. 


3.0 INTRODUCTION 


In 1986 the British Royal Aerospace Establishment (RAE), [now the Defence Research Agency (DRA)], the 
National Aeronautics and Space Administration (NASA), and the Office National D’etudes Et de Recherches 
Aerospatiales (ONERA) initiated a collaboration to compare and validate the predictive capabilities of each 
agency's computer codes that support aircraft icing technology. It was intended that this collaboration would 
yield a better understanding of the codes and suggest possible improvements to them. The three kinds of 
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codes compared were 1) water droplet trajectory calculations, 2) ice accretion modelling and 3) transient elec- 
trothermal deicer analysis. 

Each agency contributed experimental results for the code comparison exercises. The experimental data 
included droplet collection efficiencies, heat transfer coefficients, and ice accretion shapes. The data were 
obtained in the following icing facilities: the NASA Lewis Icing Research Tunnel and the NASA Heat Transfer 
Tunnel in Cleveland, Ohio, USA; the Artington Icing wind tunnel facility near Guilford in Surrey, ENGLAND; and 
the engine icing facilities at CEPr, Saclay, FRANCE. In addition to comparisons with experimental results, code 
predictions were mutually compared for hypothetical icing conditions. 

Three separate technical reports have been published, one for each of the three kinds of codes. Each 
report was authored jointly by the three agencies. DRA published and had first authorship of Part I - Droplet 
Trajectory Calculations (Ref. 1), NASA Part II - Ice Accretion Prediction (the present report), and ONERA Part 
III - Electrothermal Deicer Analysis (Ref. 2). 

These reports documented the joint studies done from 1986 to 1991, and as such, these reports provide 
the status of each agency’s code development during that time period. Since 1991, each agency has been 
learning more about their codes and continually improving them, and they have also been developing new 
codes with more capabilities (Ref. 3-10). 

This present report presents the results of the code comparison exercises for ice accretion on two dimen- 
sional airfoils. Section 4 briefly describes the main features of each agency’s ice accretion code. All three 
codes have some common features, such as, a flowfield module, a water droplet trajectory module, a control 
volume module that performs both a mass balance and a heat balance on the water and predicts ice growth. 

The flowfield solutions are obtained from well established flow codes, and the droplet trajectory calcula- 
tions are described in Reference 1 . But this is the only report that describes the heat and mass transfer calcu- 
lations. Heat transfer over ice surfaces is not well understood (Ref. 11), and therefore Section 4 presents the 
equations used to calculate the convective heat and mass transfer coefficients and the way that surface rough- 
ness is modeled. These methods vary slightly from code to code. 

Section 5 describes the experimental databases for both the heat transfer coefficients and the ice accre- 
tion shapes. Section 6 presents the heat transfer coefficient comparisons, and Section 6 also presents the ice 
shape comparisons. 


4.0 DESCRIPTION OF THE CODES 


Although each code contains the same kinds of modules, these modules may vary from code to code. For 
example, there are two different flow field modules, three different droplet trajectory integration schemes, three 
different heat and mass transfer expressions, and three different ways of predicting the ice growth during each 
time step. 

NASA 

The LEWICE code (Ref. 3) consists of three main modules as shown in Figure 1: 1) a Hess-Smith 2D 
potential flow panel method (Ref. 12) to calculate the flow field about an arbitrary two dimensional body, 2) a 
2D water droplet trajectory module (Ref. 13) that calculates the water droplet collection efficiency (flux pattern) 
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on the body, 3) a thermodynamic module that solves the mass and energy balance equations in control vol- 
umes, determines the fraction of water that freezes in each control volume, and obtains an ice profile. The 
energy balance equation, which is based on the derivation by Messinger (Ref. 14), accounts for heat and mass 
transfer by convection to the freestream, latent heat released during freezing, aerodynamic heating, and 
incoming and outgoing enthalpy of the water. Ice is assumed to grow normal to the clean airfoil or to the ice 
surface covering the airfoil. 

In LEWICE, the geometry of the 2D airfoil consists of line segments joining a discrete number of coordi- 
nates that define the surface. The length of line segments can vary around the airfoil, and smaller segments 
can be concentrated near regions where the radius of curvature is smallest. The airfoil surface is divided into 
control volumes that coincide with the line segments. 

The calculation starts at the stagnation point, where water droplets enter the control volume (c.v.) via 
impact at a rate calculated by the droplet trajectory module. A heat balance is made on the c.v., which estab- 
lishes the freezing rate of water in that c.v. Any water that does not freeze in the first c.v. is assumed to flow 
downstream to the next c.v., where water also enters by droplet impact. A heat balance is made on the second 
c.v., and any unfrozen water runs back to the third c.v. The calculation is repeated for successive control vol- 
umes until all the water that impacted the airfoil surface either freezes or reaches the trailing edge. 

Heat transfer coefficients are calculated from an integral boundary layer analysis described in Ref. 1 5. The 
laminar heat transfer coefficient does not depend upon surface roughness, but the transition point and the tur- 
bulent heat transfer coefficient are highly dependent on this variable. Transition occurs when the Reynolds 
number based on the roughness height, Re k , exceeds 600. The method used to calculate turbulent heat trans- 
fer coefficient is described in Ref. 3. 

The LEWICE code has a time-stepping capability to grow the ice. With time-stepping, ice accretes on the 
airfoil for a given time increment, and this ice growth modifies the airfoil leading edge shape. The flowfield is 
recalculated for the modified airfoil shape, and then the droplet trajectory code is rerun for the new flowfield 
and the modified airfoil shape. With the new airfoil shape and new collection efficiency as input to the code, 
another layer of ice is grown, which further modifies the leading edge. This process is repeated for enough 
time increments to reach the full icing exposure time. 

Note: Current (Ref. 10) developments have led to codes for which smaller time steps do provide improved 
ice shape definition at the cost of additional computer time. LEWICE currently uses an automated time step- 
ping procedure based on limiting the ice growth at each time step 

Inputs required by LEWICE include the airfoil geometry, cloud properties such as liquid water content 
(LWC) and droplet size (or droplet distribution), atmospheric and flight conditions such as freestream velocity 
(VJ, static temperature (T s ), static pressure (P s ), humidity (r h ), and angle of attack (a), as well as the total icing 
time. 

ONERA 


The potential equation for the flowfield is solved in a “C" grid by using a finite element method (Ref. 1 6). 
The code takes the compressibility effect into account but doesn't take the viscous effect into account. More- 
over, it cannot determine and compute a boundary layer separation. 

The trajectories are computed considering that the only forces to take into account are the drag forces, 
produced by the velocity differences between the droplet and the flowfield around it. 
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The determination of the heat transfer coefficient uses the Makkonen (Ref. 17) correlation. The roughness 
of the profile has only an effect if the boundary layer is turbulent, i.e., the Reynolds number relative to the mean 
roughness height is greater than 600. The normalized (divided by the chord, c) equivalent sand grain rough- 
ness height is a mean value always equal to 0.001 . 

The code uses the Messinger equations (Ref. 14) to determine the local freezing fraction. Starting at the 
stagnation point, this code solves the steady state heat balance equation for each successive control volume. 
It calculates the local equilibrium surface temperature, the freezing fraction and the local rate of freezing. The 
unfrozen fraction is then allowed to run back to the next control volume. 

The ONERA model can be time dependent but generally uses a different method - which is similar to a pre- 
dictor corrector method - which gives relatively good results: The ice shape is first estimated for a given icing 
time in one step.The flowfield, the trajectories and the heat transfer coefficient are then calculated on this “esti- 
mated shape". Assuming that the values of the local collection efficiency and the heat transfer coefficients are 
varying linearly from their values on the clean airfoil to their values on the profile covered by the estimated 
shape, the thermodynamic balance is made and the ice shape calculated. The codes are used twice for all the 
icing durations. 

Inputs required by ONERA include the airfoil geometry, cloud properties such as liquid water content 
(LWC) and droplet size (or droplet distribution), atmospheric and flight conditions such as freestream velocity 
(V„,), static temperature (T s ), static pressure (P s ), and angle of attack (a), as well as the total icing time. 

DRA 

The DRA ice accretion code, called TRAJICE2, is structured very similar to LEWICE. TRAJICE2 (Ref. 4) 
uses a two dimensional flow code developed by Fiddes (Ref. 18) coupled with a 2D particle trajectory code to 
determine the mass of water impacting on the surface. A thermodynamic module evaluates the governing heat 
and mass balance equations of an iced surface to determine the rate of ice accretion from which an ice profile 
can be obtained. The heat balance adopted in the DRA code is based on the standard Messinger approach but 
with both the convective cooling and the evaporative cooling terms modified to allow for compressibility effects. 
In particular, the changes make the evaporative cooling term more sensitive to the local surface pressure than 
in the standard Messinger expression. The airfoil is divided into equal-spaced elements for which the ice 
growth is assumed to be normal to the surface for glaze ice or parallel to the freestream for rime ice. TRAJICE2 
is also a time stepping code capable of predicting the ice shape using a single or multiple step process. Heat 
transfer coefficients are calculated using an integral boundary layer routine similar to LEWICE if it is a multiple 
time step run or by using correlations for cylinders when using a single time step where the cylinder diameter is 
an input variable and not necessarily equal to the airfoil leading edge diameter. TRAJICE2 also uses a 50% 
higher heat transfer coefficient when comparing to tunnel test to account for the higher turbulence in wind tun- 
nels. 

Although using multiple time steps should, in principle, be more accurate than a single time step, it was 
found by DRA that single time steps sometimes gave better comparisons with experiment than did multiple 
time steps. Furthermore, no determination could be made if there was an optimum time step increment (such 
as 2 min, 1 min, 30 sec, 20 sec, etc). 

Inputs required by TRAJICE2 include the airfoil geometry, cloud properties such as liquid water content 
(LWC) and droplet size (or droplet distribution), atmospheric and flight conditions such as freestream velocity 
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(VJ. static temperature (T s ), static pressure (P s ), and angle of attack (a), as well as the total icing time (Ref. 
4,5). 

Heat Transfer Analysis 

Before the heat transfer comparisons are presented, a brief overview of the heat transfer analysis is pre- 
sented below. As mentioned previously, NASA and ONERA use the integral boundary layer method exclusively 
to calculate h, while DRA uses correlations for a cylinder for single time step runs and the integral boundary 
layer method for multiple time step runs. However, there are a few differences among them which will be noted. 
It is not the intent of this paper to show the entire derivation of the heat transfer analysis, but merely to highlight 
the significant terms that affect the calculation of h and to show how the turbulent heat transfer coefficient is 
highly dependent on the local roughness. A complete derivation of the heat transfer analysis is contained in 
Appendix B of Ref. 3. 

Each code uses a control volume approach to develop its mass and energy balance equations. An airfoil or 
other body is defined by a discrete number of points which are connected by line segments, each segment 
defining a control volume. Figure 2 illustrates the mass balance for a clean or iced surface which accounts for 
water entering and leaving a control volume through impingement of the particles, evaporation, and water from 
one control volume to the next. The energy balance terms, shown in Figure 3, account for conduction and con- 
vection to the freestream, evaporation and aerodynamic heating. The surface roughness, k, is expressed as an 
equivalent sand grain roughness (more will be said on this later). A thorough discussion of the control volume 
concept used to formulate the mass and energy balance equations is found in Appendix A of Ref. 3. 

The first step in the heat transfer analysis is to define the transition location from laminar to turbulent flow. 
This is computed by using the following empirical relationship for the local roughness Reynolds number 


Re k = 



(i) 


where V k is the velocity at the roughness height, k is the equivalent sand grain roughness and u is the vis- 
cosity. For laminar flow, Re k will be less than 600. When Re k is greater than 600 the flow is assumed to be tur- 
bulent. 


Once the transition location is determined, the heat transfer coefficient can be determined for each seg- 
ment of the airfoil. For laminar flow, Re k less than 600, the local heat transfer coefficient is calculated using the 
following equation developed by Smith and Spaulding (ref. 19) 


, 0.296 X. 

h, = — — 

cjRe 




(rO 


n 1/2 


- 0 


( 2 ) 


where X is the thermal conductivity and V e is the velocity at the edge of the boundary layer calculated from 
the potential flow code. 
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As shown, the laminar heat transfer coefficient is not dependent on the local roughness. However, for tur- 
bulent flow, the heat transfer coefficient is directly dependent on the local roughness and is calculated using a 
technique taken from Ref. 1 5 for turbulent flow 


Nu t = 


\ c f Re x Pr 


( 3 ) 


where Re x is the local Reynolds number, St k is the roughness Stanton number, Pr is the Prandtl number, 
Pr t is the turbulent Prantl number, and c f is the skin friction coefficient. As will be shown later, the value of the 
roughness, k, can significantly influence h. 


To calculate the turbulent heat transfer coefficient one must first calculate the skin friction coefficient and 
the roughness Stanton number which are also dependent on the local roughness. 


The skin friction coefficient is determined by the following equation 



0.41 


, ( 8649 

l ,0 <— + 2 


.586) j 


(4) 


where 0, is the turbulent momentum thickness which can be evaluated by Thwaites formula (Ref. 15). 
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( 5 ) 


The roughness Stanton number is determined by the following equation (Ref. 19) 


St k = 1.92 Re^Pr** ( 6 ) 

The equivalent sand grain roughness is the only variable needed to calculate the heat transfer coefficient 
that cannot be calculated readily from boundary layer theory. Because of the irregular surface roughness ele- 
ments associated with an ice surface, none of the existing empirical correlations for equivalent sand grain 
roughness are valid. Therefore, there have been several correlations which have been developed which are 
specific to roughness elements on an ice surface. These are described below. 

A correlation relating the sand grain roughness to several icing parameters was developed by Ruff and can 
be found in Ref. 3, Appendix F. This roughness correlation is currently used in LEWICE and was developed in 
the following manner. To begin with, a set of ice shapes were selected from data taken on a NACA 0012 airfoil 
by DRA in the Artington icing wind tunnel facility. Next, LEWICE was used to predict the ice shapes for each 
case by arbitrarily varying the roughness until the best predicted ice shape was obtained. A correlation that 
relates the equivalent sand grain roughness to the freestream velocity V„, liquid water content (LWC) and the 
static temperature (T s ) was then developed. The final form of the correlation for the equivalent sand grain 
roughness is: 
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C V. C Jbase\ C J LWC\ C/T\ C/V 


( 7 ) 


where c is the airfoil chord and (k/c) base is a baseline value accounting for characteristics that are unique to 
the facility in which the ice shapes were formed. 


Note: The current (Ref. 10) version of LEWICE predicts roughness height from the surface physics, and is 
not based on the correlation above. See Ref. 8 for a full description. The equation used is 


k = 


(a 2 W 3 

4op 

KpF X > 


( 8 ) 


where a is the surface tension, p. is the viscosity of water, p is the water density, x is the surface shear 
stress and F is the fraction of the surface which is wetted. Reference 10 shows good comparisons of this 
model to measured ice roughness heights. 


The DRA uses one of two methods to establish a value for the surface roughness parameter. If an equiva- 
lent sand grain roughness can't be estimated, the maximum roughness height of the elements can be used. 
Typical values range from 0.1 to 5 mm and are input by the user. If a negative value of roughness is specified, 
the code will automatically calculate the equivalent sand grain roughness by using an empirical relationship 
similar to that used by LEWICE which is also a function of V„,, LWC, and T s . TRAJICE2 also uses a turbulence 
level correction factor to calculate the heat transfer coefficient. A multiplication factor of 1.0 or 1.5 is recom- 
mended depending on whether the comparisons are made to natural or artificial (i.e., tunnel) icing data, 
respectively (refs. 4,5). This correction factor is attributed to the higher turbulence levels in icing tunnels. In 
addition, the DRA uses the cylinder correlation (which requires a single time step) for most of the glaze ice 
shape predictions shown in this report and uses the integral boundary layer approach (which can use multiple 
time steps) for the rime cases and the warmest glaze icing conditions. 

Initially, k is used to calculate the turbulent Reynolds number for each segment to determine the point of 
transition from laminar to turbulent flow as governed by Equation 1. If the turbulent Reynolds number is less 
than 600, the heat transfer coefficient is calculated using Equation 2 which is a function of the edge velocity 
only. If Re k is greater than or equal to 600, the heat transfer coefficient must then be calculated using Equation 
3 which is highly dependent on the roughness. 


5.0 EXPERIMENTAL DATABASE 


Below is a brief description of both the heat transfer and ice shape database selected for the DRA/NASA/ 
ONERA collaboration on icing research. 

Heat transfer data was taken on four simulated ice shapes which were based on shapes measured in the 
NASA Lewis 6' X 9' Icing Research Tunnel (IRT) as shown in Figure 4 (ref. 11). The simulated ice shapes were 
produced by numerically smoothing the original ice shape tracings taken on 2" cylinders to create a symmetri- 
cal ice shape about the cylinder axis. The simulated ice shapes were then fabricated out of polyurethane foam 
and fitted with electrically heated copper strips which were embedded in the foam at closely spaced intervals 
around the surface of the ice shape. Each copper strip contained a thermocouple which was used to ensure 
the temperature of the surface was isothermal. The amount of power required to maintain a given surface tem- 
perature was then used to determine the heat transfer coefficient of the surface. 
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Over the years, NASA has compiled a database of ice shapes on various airfoils and other configurations 
(refs. 20, 21, 22) which have been used for code validation purposes. The majority of ice shape data collected 
by NASA has been on 21” chord NACA 0012 airfoils. Since temperature is a very important parameter that 
affects the type and shape of an ice accretion, the ice shape data in Figure 5 (ref. 20) was selected to examine 
the performance of the three codes as a function of temperature. The two sets of ice shape data, shown in Fig- 
ure 5, were taken on a NACA 0012 airfoil at 4° angle of attack with the indicated icing conditions. This figure is 
reproduced from an earlier report (Ref. 21). Not all of the ice shapes shown in this figure were included for this 

comparison. The first set of ice shapes are for a velocity of 58 m/s and an LWC of 1 .30 g/m 3 . The velocity was 
increased to 94 m/s with an LWC of 1 .05 g/m 3 for the second set of ice shapes. As shown in Figure 5, the effect 
of temperature on the ice shape is clearly seen as the total temperature is increased from -26°C to 0°C. At the 
coldest temperatures, the ice shapes are completely rime and are streamlined similar to the clean airfoil. The 
ice shapes become a mixture of both rime and glaze ice with the glaze ice horns becoming more prevalent as 
the temperature is increased. For temperatures between -26°C and - 5°C the mass of ice remains relatively 
constant. For temperatures warmer than -5°C the total mass of ice accreted begins to decrease until, at 0°C, 
there is no ice accreted. The decrease in the mass of ice accreted is attributed to water blow-off or water shed- 
ding at the warmer temperatures. It is this temperature region that ice prediction codes have difficulty since the 
modeling of water shedding is either neglected or only approximated in the present models. Note: The current 
(Ref. 10) version of LEWICE has this capability. 

An SA13112 airfoil was used by ONERA to obtain experimental ice shapes for various Mach numbers and 
angles of attack (Ref. 23). From these data, several cases were selected and are shown in Figure 6a. For a 
Mach number of 0.8, code comparisons were made even though no experimental data was reported by 
ONERA. 

Additional ice shape data for a NACA 0012 airfoil has been reported by the DRA (refs. 24, 25). From these 
data, several experimental ice shapes were selected for comparison purposes and are shown in Figure 6b. Ice 
shapes from all three icing facilities were included in this study to determine the effects, if any, of comparing the 
codes to ice shape data taken from differing icing tunnel facilities. 


6.0 RESULTS AND DISCUSSION 


The results will be separated into two sections. The first section will present a brief overview of the integral 
boundary layer method used by all three codes to calculate the heat transfer coefficient (all three ice accretion 
codes use essentially the same equations to calculate h although there are some differences which will be 
noted). The calculated heat transfer coefficients from the three ice accretion codes will be compared to experi- 
mental values obtained for the four simulated ice shapes shown in Figure 7. The second section will compare 
ice shape predictions from each of the three ice accretion codes to the experimental ice shape data shown in 
Figures 5 and 6. 

The four simulated ice shapes shown in Figure 7 were used to obtain heat transfer data and were reported 
by Van Fossen, et al. (Ref. 11). Heat flux gauges were located at the positions shown for each shape. For each 
of the simulated ice shapes, heat transfer data were taken for both a smooth and roughened surface. A highly 
polished surface was used for the smooth surface while the rough surface was simulated by introducing grains 
of sand onto the leading edge of the simulated ice shape surface. 

Figures 8 and 9 show the comparison of the heat transfer coefficients calculated by each of the three 
codes for the 5 minute simulated glaze ice shape with both a smooth and rough surface, respectively. As 
shown, all three codes predict the heat transfer coefficient reasonably well in the stagnation region, but signifi- 
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cantly overpredict h near the horns of the ice shape. Aft of the horn region, all three codes do not compare well, 
probably due in part to inaccurate modeling of the flowfield. Aft of the hom region, there is a large separation/ 
reattachment bubble which has been examined through calculations made with a Navier-Stokes code called 
ARC2D (Ref. 6). Since none of the codes model the flow separation/reattachment, the resulting calculated flow 
velocity is larger than the actual velocity in this region. 


Figures 10, 11 and 12 show the calculated pressures, velocities and temperatures at the edge of the 
boundary layer versus the non-dimensional surface distance from the stagnation point for the 5 minute glaze 
ice shape. As shown, all three codes calculate very similar boundary layer properties. It is interesting to note 
that each code predicts a large increase in the velocity near the horns of the ice shape while the pressure and 
temperature plots indicate a significant decrease in pressure and temperature near the horns. This is due to 
the large drop in the calculated pressure coefficient (not shown) in the horn region which can be attributed to 
the inaccurate modeling of the flow separation that occurs. Therefore, the calculated flow remains attached 
and must increase significantly to turn past this region. The net result is significantly higher heat transfer being 
calculated in this region. 


Figure 13 shows the LEWICE calculated heat transfer coefficient for the 5 minute glaze ice shape with a 
rough surface (k=0.33 mm). For s/c values between 0.0 and 0.53, Re k is less than 600 and the heat transfer 
coefficient is calculated using Equation 2. In this region, the calculated h compares reasonably well with exper- 
iment For s/c values between 0.53 and 2.10, Re k is greater than 600 and h is calculated using Equation 3. 
Notice the large spike in the calculated h at s/c=0.53. This is the point where the calculated flow transitions 
from laminar to turbulent and is evident by the large increase in h from around 100 to 1050 W/m 2o K. For s/c 
values between 0.53 and 2.10, the calculated h does not compare well with experiment The same calculation 
was repeated with LEWICE using a k=1 .5*1 O' 3 mm to simulate a smooth surface. In this case, Re k is less than 
600 over the entire ice shape and h is calculated using Equation 2. Therefore, by merely increasing or 
decreasing the value of k, one can significantly affect the point of transition from laminar to turbulent as well as 
the magnitude of h. 

Figures 14 and 15 show the experimental and calculated heat transfer coefficients for the 15 minute simu- 
lated glaze ice shape with a smooth and rough surface, respectively. In both instances, the calculated h's com- 
pare reasonably well in the stagnation region. In the horn region, DRA and ONERA codes moderately 
overpredict h, while the NASA calculation which underpredicts h for the smooth surface. Figures 16, 17 and 18 
are the calculated values of P e , V e and T e for the 1 5 minute glaze ice shape. As shown, the three codes predict 
very similar boundary layer properties for P e and V 0 . The NASA results forT e are much higher than the results 
reported by DRA and ONERA. In the horn region of the ice shape, there is a large increase in the velocity and 
a large decrease in the pressure and temperature as was observed previously for the 5 minute glaze ice 
shape. 

The calculated transition point for the 15 minute glaze ice shape with a rough surface (k=0.33 mm) is 
shown in Figure 1 9. In this case, the calculated transition point is at s/c=0.927. Once again, NASA predicts h 
reasonably well up to the point of transition. Aft of this region, the calculated h does not compare well with 
experiment. The calculations were repeated using a k=1.5*10' 3 mm but no transition point was calculated. 
Once again, the calculated flow remained laminar over the entire ice surface by using a small value of k. 

Heat transfer coefficients for the 15 minute simulated rime ice shape are shown in Figs. 20 and 21 for both 
a smooth and rough surface, respectively. For the smooth surface, each code significantly underpredicts h up 
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to the point of flow separation which occurs at the top and bottom of the cylinder. For the rough surface, each 
code predicts h very well along the entire iced surface up to the point of separation. Figures 22, 23 and 24 
show the P e , V e and T e calculated by each code for the 15 minute rime ice shape. Once again, each of the 
three codes predict very similar boundary layer properties. The calculated transition point for the 15 minute 
rime ice shape is shown in Fig. 25 for a rough surface. In this case, the transition point occurs at s/c=0.161 
which is much closer to the stagnation point than the other two cases previously examined. However, unlike 
the other two cases, LEWICE predicts h in both the laminar and turbulent regions well. 

As mentioned previously, DRA uses the cylinder values for many of the glaze ice predictions while NASA 
and ONERA use the integral boundary layer method throughout. DRA shows that the ice shape prediction is 
better when the cylinder correlation (a user-selected procedure) is used. Figure 26 shows a comparison of the 
heat transfer coefficient produced on a NACA0012 airfoil using the integral boundary layer technique and the 
heat transfer coefficient on a cylinder . The cylinder diameter is chosen so that it is the same as the leading 
edge diameter of a NACA0012. As seen in this figure, the heat transfer coefficient is much higher for the cylin- 
der. This is to be expected, since cylinder values are symmetric whereas the boundary layer technique gives 
lower values due to the lower velocities on the pressure surface of the airfoil. Although the courseness of the 
experimental data makes any final conclusion open to interpretation, the previous examples show that the inte- 
gral boundary layer method either predicts or over predicts heat transfer coefficients regardless of which code 
is used. The cylinder results as shown in Figure 26 exceed the IBL results, so they would be in the direction 
opposite what is needed to predict the correct HTC value. It appears to illustrate that there is some term miss- 
ing from all of the codes which would account for this additional heat loss. 

Ice Shape Prediction 

Each of the three ice accretion codes were used to predict the ice shapes for a selected matrix of experi- 
mental ice shapes taken on both a NACA 0012 and an SA13112 airfoil. For each of the 16 cases examined, 
the experimental ice shapes are indicated with a dashed line and the predicted ice shapes are indicated by a 
solid line. Table 1 shows a summary of all of the conditions for these cases 

Figures 27-36 show comparisons between the experimental and predicted ice shapes for the experimental 
data of Olsen (Ref. 20) previously shown in Figure 5. In each case, the predicted ice shapes show excellent 
comparison to experiment except for the mixed glaze condition in Fig. 29. 

Figures 27 and 28 show the predicted ice shapes for a velocity of 58 m/s and total temperatures of -26° 
and -18°C, respectively. These shapes are rime ice shapes, although the ONERA code shows the onset of a 
small horn formation not seen on the experiment shape or the shape predicted by DRA and NASA. 

Figure 29 shows the comparison at -12°C which is a mixed glaze condition. A mixed glaze condition is a 
transitional case from rime to glaze where the center is glaze ice, but the horns are large rime feathers. In this 
case, the predicted ice shapes compare reasonably well in terms of ice mass but tend to miss some features of 
the experimental ice shape. DRA predicts a rime shape for this condition, while NASA and ONERA predict a 
classic double horn glaze shape. This difference in ice type shows the typical difficulty in predicting this class of 
shapes. 

Figure 30 shows the predicted ice shapes for a total temperatures of - 5°C. In this case, the predicted ice 
shapes are similar to each other and compare well with the experimental ice shapes. DRA shows a larger 
lower horn than the experiment or the other two codes. It should be noted that DRA used a single time step 
and a cylinder correlation for this case, which produces a much higher value of heat transfer coefficient on the 
lower surface as was shown in Figure 26. 
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Figure 31 shows the prediction at -3°C total temperature. Here, the predictions are smoother than the 
experimental shape. The predicted horns are not as well developed and the upper impingement is slightly 
overpredicted, but still the overall shape of the prediction is still very good. 

Figure 32 shows the prediction at -1°C total temperature. In this case, each code overpredicts the upper 
impingement limit even though the general shape and the mass of ice agree reasonably well except for the 
additional runback ice on the upper surface. This case may illustrate the need to include improved water run- 
back models which include shedding of excess water from the surface. 

The second set of ice shape data, shown in Figure 5, is for a velocity of 94 m/s. Predicted ice shapes for a 
total temperature of -26°C are shown in Figure 33. DRA compares reasonably well with the experimental ice 
shape except that they do not predict the pointed leading edge seen on the experimental ice shape. ONERA 
underpredicts the mass of ice for this case while NASA overpredicts the lower impingement limit and underpre- 
dicts the stagnation thickness of ice. For both ONERA and NASA, this is poor prediction for a rime ice case. 

Ice shape predictions for a total temperature of -12°C are shown in Figure 34. As opposed to Figure 30, 
where the cylinder correlation over predicted the horn size, in this case DRA obtained excellent results using 
its cylinder correlation for heat transfer instead of the integral boundary method. NASA does not predict the 
correct size of the lower ice horns although it captures the feathers aft of the horn well, while ONERA does not 
predict the hom angles well. Whether one should use the cylinder correlation or the integral boundary layer 
approach depends upon the experience of the user, and is not a true reflection of a difference in the codes. 

The logic DRA recommends for choice of equivalent cylinder or integral boundary layer method for calcu- 
lating the heat transfer coefficient distribution is in the first instance to use the cylinder method and single time 
step approach. If, based on the experience of the user, the predicted shape is unrealistic, then further calcula- 
tions with the IBL method and multi-stepping approach may be required. There is however at this time no sim- 
ple procedure as to which method should be used. 

Figures 35 shows the predicted ice shapes for total temperatures of -8°C. All three codes predict the upper 
hom accurately. DRA predicts the lower horn as well where both NASA and ONERA underpredict the lower 
hom size. Again, this is due to the use of the cylinder correlation and hence is a function of the user’s knowl- 
edge and experience and is not a reflection of the code capabilities. The cylinder correlation produces much 
higher heat transfer coefficients on the lower surface as shown in Figure 26. However as shown in Figs. 8,9, 
and 13-15, the integral boundary layer method overpredicts experimental heat transfer results. This discrep- 
ancy illustrates some physical phenomena which is not properly handled by any of the codes. 

Figure 36A shows the prediction at -2°C total temperature. The upper hom angle is well predicted by the 
codes, but there are some small differences in the size of the hom as the ONERA prediction gives slightly 
larger horns while the DRA predicts slightly smaller horns than the experiment. Once again, DRA shows supe- 
rior prediction on the lower hom due to the higher heat transfer coefficient produced by the cylinder correlation. 
Both NASA and ONERA predict a smaller hom with more runback refreezing. Once again, since an iced airfoil 
does not resemble a clean cylinder, this result shows that there is some physical phenomena which is not 
properly handled by any of the codes. Figure 36B shows a more direct comparison of the codes, as DRA used 
its integral boundary layer method instead of the cylinder correlation. In this case, the DRA prediction shows 
the same characteristics as the NASA and ONERA predictions. 

Figures 37-39 show comparisons between the predicted and experimental ice shapes for the DRA data 
(shown in Figure 6) taken on a 0.152 m chord NACA 0012 at 0°C and 8.5° angle of attack and for a beak ice 
shape on a 0.465m chord NACA 0012 at 8° angle of attack. As shown in Figure 37, for the 0° angle of attack 
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case, all of the codes underpredict the leading edge thickness of ice shown in the experiment, with NASA 
underpredicting the shape by the greatest amount. Due to the large, thin rime feathers on the experimental ice 
shape, this may be due to the ice density values used by the codes. Rime ice can have ice densities as low as 
100 kg/m 3 as opposed to a value of 917 kg/m 3 for glaze ice. DRA again increases the heat transfer coefficient 
by using its cylinder correlation. NASA and ONERA predict similar thicknesses at the stagnation point which is 
approximately half of the ice growth shown for the experimental ice shape. ONERA predicts much larger horns 
than the other codes or the experiment. The predicted horn angle is larger as well. NASA underpredicts the 
entire mass. Since ONERA is modelling more of the compressible nature of the flow at this high mach number 
case, it should produce an ice shape closer to the experiment. 

Figure 38A shows the comparison for the 8.5° angle of attack case. On the uppdr surface, all three codes 
predict a lower horn angle than the experiment and a further extent of ice impingement than the experiment. 
Due to the high angle of attack, DRA suggests that the flow was stalled on the airfoil for this case. Stall behav- 
ior can not be modeled by any of the flow codes, but its effect can be approximated by running the case at a 
lower angle of attack. DRA ran the same conditions in Figure 38A, but changed the angle of attack until the 
predicted ice shape matched the ice shape in the experiment. This result is shown in Figure 38B where a 2° 
angle of attack was used. 

The improved agreement shown in Figure 38B between the predicted results and experiment supports the 
DRA theory that the airfoil may have stalled. At present with the existing codes, there is no way of predetermin- 
ing if the ice accretion will cause the airfoil to stall. Hence all users of the codes should be aware of the poten- 
tial danger of significant errors in the predicted ice shapes when making a prediction on an airfoil at a high 
angle of attack. Additionally, DRA used a single time step and a cylinder correlation for heat transfer coefficient 
to predict the shape shown in Figure 38A, but used multiple time steps and the integral boundary layer method 
for the shape in 38B as these produced the best prediction each time. Again, it would be advantageous if a 
method existed which could ascertain which model to use for design conditions for which there is no experi- 
mental data. 

Figure 39 shows the comparison for a beak ice condition at 8° angle of attack and a total temperature of 
4.7°C. There is no experimental data at this condition. All of the codes predict similar ice shapes, with ONERA 
predicting the largest beak horn and DRA predicting the furthest aft extent. 

Figure 40 shows the predicted ice shapes for a 0.6 m chord SA13112 airfoil at M=0.25 and 10° angle of 
attack. As shown all three codes do a poor job of predicting the ice shape. However, there are several similari- 
ties between the three predicted ice shapes. First, each code predicts a large horn near the leading edge of the 
airfoil. Second, each code predicts a fair amount of ice on the lower surface of the airfoil. The predicted ice 
shapes are typical ice shapes that would be expected for an angle of attack of 10° yet the experimental ice 
shape looks more like an ice shape accreted at 2°. This would again indicate a condition for which the flow has 
stalled. DRA again parameterized angle of attack and again determined that the angle of attack which best fits 
the experiment is 2°. This shape is shown in Figure 40B. DRA was consistent in again using a single time step 
for the original angle of attack and multiple time steps for the stalled angle of attack. Note that although the 
case shown in Figure 38A was an 8° angle of attack condition and the case shown in Figure 40A is a 10° angle 
of attack condition, the stalled angle of attack was again 2°. Higher accuracy flow codes such as Navi6r-Stokes 
would be needed to better predict these ice shapes without artificially reducing the angle of attack. 

Figure 41 shows the predicted ice shapes for a SA13112 airfoil at M=0.50 and 0° angle of attack. All three 
codes predict ice much further back than shown on the experiment. This can again be attributed to the use of 
low complexity flow solvers which allow particles to impinge aft of ice horns. Both NASA and ONERA predict a 
clean leading edge while DRA matches the stagnation ice thickness due to their use of higher heat transfer 
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coefficient values. Both DRA and ONERA predict the horn size, with a horn angle which is slightly larger than 
the experiment while NASA underpredicts the ice shape throughout. 


Figure 42 shows the predicted ice shapes for a SA13112 airfoil at M=0.8 and 10° angle of attack for which 
there is no experimental data. For this case, all three codes show similar double horn glaze ice shapes. Simi- 
larly to the other conditions, NASA and ONERA match closely in the stagnation region, while DRA has more 
ice throughout 

7.0 CONCLUSIONS/RECOMMENDATIONS 


Each of the three ice accretion codes were used to predict ice shapes for a selected matrix of experimental 
ice shapes. All three ice accretion codes did a reasonable job of predicting the ice shapes. For each of the 
experimental ice shapes examined, one of the three ice accretion codes (LEWICE, ONERA or TRAJICE2) pre- 
dicted the general features of the experimental ice shape (i.e. shape, impingement limits, mass of ice) better 
than the other two. Additional comparisons were made to heat transfer data taken on three simulated ice 
shapes. All three codes were judged to compare reasonably well with the experimental heat transfer data in the 
stagnation region of the ice shapes. In the horn region of the ice shapes, all three codes tend to overpredict the 
heat transfer coefficient, sometimes quite significantly. Several reasons for this may be the increased velocities 
calculated by the potential flow code around the horn region and/or an inaccurate roughness correlation con- 
tained in the model. 

The results of the codes can also be improved significantly by altering the user input, as indicated by the 
difference between the one step and multi-step results provided by the DRA. For this reason, this code 
requires a degree of user experience in order to achieve the optimum accuracy. Also, the codes used in this 
report are unlikely to provide good agreement with experimental ice shapes formed on an airfoil at a high 
angle of attack (e.g. £ 8°) because aerodynamic stall caused by the ice accretion is neglected in these models. 
Selective use of cylinder correlations despite concerns about the accuracy of the resulting heat transfer coeffi- 
cient, as well as use of reduced angles of attack in stall conditions can greatly improve the predictions of the 
ice shape. Unfortunately, there is no methodology for determining the use of either correction for cases where 
the ice shape is not known beforehand. 

The following recommendations are made in support of an ongoing collaboration between NASA/ONERA/ 
DRA to validate an upgrade future ice accretion codes: 

1 . Additional comparisons be made to ice shape data taken on other airfoils and from other icing facilities. 

2. Additional experimental heat transfer data is needed to assess the accuracy of the existing codes. 

3. Quantitative data on ice roughness is needed to establish better estimation of the equivalent sand grain 
roughness parameter k. 

4. Additional research into the aerodynamic effects of icing which can then be used for future improvement of 
ice accretion models. 
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No Experimental Data Available 

















Figure 4 - NASA Lewis Icing Research Tunnel 


O AIRSPEED. 209 km/hr; LWC, Ug/m 3 ; TIME, I min. 


TOTAL 

TEMPERATURE 


□ 



-26° C -20° C -18° C -15° C -12° C -8° C -5°C -2°C -1°C 0°C 

AIRSPEED. 338 km/hr; LWC. L05g/m 3 ; TIME, 6.2mln. 


TOTAL 

TEMPERATURE 



Figure 5 - Summary of NASA Ice Shapes 
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Figure 6(a) ONERA Icing Data Figure 6(b) DRA Icing Data 



(c) Fifteen minute glaze ice (d) Fifteen minute rime ice 


Figure 7 - Ice shape cross sections showing profile coordinates and 

heat fluz gage locations 
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Figure 8 * 5 min Glaze Iced Cylinder 
(Smooth Case) 
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Figure 9 • 5 min Glaze Iced Cylinder 
(Rough Case) 
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figure 12 - 5 mb Glaze Iced Cylinder 
(Smooth Case) 
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Figure 15 - 15 min daze Iced Cylinder 
(Rough Case) 
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Figure 18 • 15 min daze Iced Cylinder 
(Smooth Case) 
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Figure 19 • 15 min daze Ice 
(Rough Case) 
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Figure 20 - 15 min Rime Iced Cylinder 
(Smooth Case) 
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Figure 21 • 15 min Glaze Iced Cylinder 
(Rough Case) 
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Figure 22 • 15 min Rime Iced Cylinder 
(Smooth Case) 
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Figure 23-15 min Rime Iced Cylinder 
(Smooth Case) 
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Figure 24 - 15 min Rime Iced Cylinder 
(Smooth Case) 
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Figure 25 - 15 min Rime Ice 
(Rough Case) 
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Figure 26 - Comparison of Heat Transfer 
Coefficient for a NACA0012 at a = 4* 
and a Cylinder 
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Figure 27 - Comparison of ice shape predictions for NASA cases. 

NACA0012, Chord = 0.53, a = 4.0, V® = 58.1, Ts = -27.8‘C, LWC = 1.3, MVD = 20.0 M , r = 480.0 
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Figure 28 - Comparison of ice shape predictions for NASA cases. 

NACA0012, Chord = 0.53, a = 4.0, V® = 58.1, Ts = -19.8‘C, LWC = 1.3, MVD = 20.0fi, r = 480.0 
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Figure 29 - Comparison of ice shape predictions for NASA cases. 

NACA0012, Chord = 0.53, a = 4.0, V® = 58.1, Ts « -13.9*C, LWC = 1.3, MVD = 20.0^, r = 480.0 
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Figure 30 - Comparison of ice shape predictions for NASA cases. 

NACA0012, Chord = 0.53, <r = 4.0, V® = 58.1, Ts = -6.7°C, LWC = 1.3, MVD = 20.0/x, r = 480.0 
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Figure 31 - Comparison of ice shape predictions for NASA cases. 

NACA0012, Chord = 0.53, a = 4.0, V« = 58.1, Ts * -3.9*C, LWC = 1.3, MVD = 20.0^., r = 480.0 


36 




Experiment 

Prediction 



DRA 

1 Time Step 



NASA 

4 Time Steps 



Figure 32 - Comparison of ice shape predictions for NASA cases. 

NACA0012, Chord = 0.53, « = 4.0, V« = 58.1, Ts = -2.8*C, LWC = 1.3, MVD = 20.0/x, r = 480.0 
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Figure 33 - Comparison of ice shape predictions for NASA cases. 

NACA0012, Chord = 0.53, a = 4.0, V® = 93.9, Ts « -30.5*C, LWC = 1.05, MVD = 20.0/t, r = 372.0 
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Figure 34 - Comparison of ice shape predictions for NASA cases. 

NACA0012, Chord = 0.53, « = 4.0, V® = 93.9, Ts = -16.6 0 C, LWC = 1.05, MVD = 20.0^, r = 372.0 
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Figure 35 - Comparison of ice shape predictions for NASA cases. 

NACA0012, Chord = 0.53, or = 4.0, V® = 93.9, Ts = -12.2‘C, LWC = 1.05, MVD = 20.0^, r = 372.0 
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Figure 36A - Comparison of ice shape predictions for NASA cases. 

NACA0012, Chord = 0.53, a = 4.0, V® = 93.9, Ts - -6.6‘C, LWC = 1.05, MVD = 20.0^, r = 480.0 
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Figure 36B - Comparison 

NACA0012, Chord = 0.53, a « 4.0, V 


of ice shape predictions for NASA cases. 

® = 93.9, Ts = -6.6‘C, LWC = 1.05 MVD » 20.0/x, r = 480.0 
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Figure 37 - Comparison of ice shape p: 
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jdictions for DRA cases. 

LWC = 0.5, MVD = 20.0/i, r = 480.0 
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Figure 38B - Comparison of ice shape predictions for DRA cases. 

NACA0012, Chord = 0.53, a = 2.0, V» = 130.5, Ts = -12.3°C, LWC = 0.5, MVD = 17.5/1, r = 120.0 
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Figure 39 - Comparison of ice shape predictions for DRA cases. 

NACA0012, Chord = 0.53, a - 8.0, V® = 131 . 5 , Ts = -3.9‘C, LWC « 0.6, MVD = 20 . 0 / 1 , r = 180.0 
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Figure 40B - Comparison of ice shape predictions for ONERA cases. 

SA13112, Chord = 0.6, a = 2.0 (stalled), M = 0.25, Ts = -10.0°C, LWC = 0.5, MVD = 20.0 ^ r = 900.0 
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Figure 41 - Comparison of ice shape predictions for ONERA cases. 

SA13112, Chord = 0.6, a = 0.0, M = 0.50, Ts = -10.0‘C, LWC = 0.5, MVD = 20.0^, r = 450.0 
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Figure 42 - Comparison 

SA13112, Chord « 0.6, a ■* 0.0, 


of ice shape predictions for ONERA cases. 

M = 0.80, Ts = -30.0'C, LWC = 0.5, MVD = 20.0 M , r = 180.0 
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